library(lme4)
library(lmerTest)
library(dada2)
library("phyloseq")
library(ggplot2)
library(tidyverse)
library(fantaxtic)
library(DECIPHER)
library(phangorn)
library(reshape)
library(janitor)
library(multcomp)
library(MuMIn)
library(SciViews)
library(mblm)
library(regclass)
library(devtools)
library(broom)

# define predictive R^2 functions for later
PRESS <- function(linear.model) {
  #' calculate the predictive residuals
  pr <- residuals(linear.model)/(1-lm.influence(linear.model)$hat)
  #' calculate the PRESS
  PRESS <- sum(pr^2)
  
  return(PRESS)
}

pred_r_squared <- function(linear.model) {
  #' Use anova() to get the sum of squares for the linear model
  lm.anova <- anova(linear.model)
  #' Calculate the total sum of squares
  tss <- sum(lm.anova$'Sum Sq')
  # Calculate the predictive R^2
  pred.r.squared <- 1-PRESS(linear.model)/(tss)
  
  return(pred.r.squared)
}

# set scientific notation off
options(scipen=999)

# load newest data
load("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code/functional type by sample.RData")

# load phyloseq object and exploration type and guild data
#load("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code/RA otu by guild and exploration type.RData")


#######################
# SAPROTROPHS #
############################

# relative abundance
all.parms <- lm(saprotroph ~ organic.soil.C.N + mineral.soil.C.N + organic.soil.GWC + slope + ln(Distance) +
                  carcass.load + 
                  ph + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(formula = saprotroph ~ organic.soil.C.N + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# saprotroph relative abundance decreased with organic soil C:N

# richness
all.parms <- lm(saprotroph.richness ~ organic.soil.C.N + mineral.soil.C.N + GWC + slope + ln(Distance) + ph + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(mineral.soil.C.N=mean(exploration.guild.RA$mineral.soil.C.N), 
                     GWC=mean(exploration.guild.RA$GWC), 
                     decomposing.carcass="no")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")


step.model.lmer <- lmer(formula = saprotroph.richness ~ mineral.soil.C.N + decomposing.carcass + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# saprotrophic richness increases with decomposing carcass and decreases with mineral soil C:N

# diversity
exploration.guild.RA.sap.div <- exploration.guild.RA %>% drop_na(saprotroph.diversity)
all.parms <- lm(saprotroph.diversity ~ organic.soil.C.N + mineral.soil.C.N + GWC + slope + ln(Distance) + carcass.load + ph + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA.sap.div)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(organic.soil.C.N=mean(exploration.guild.RA$organic.soil.C.N),
                     mineral.soil.C.N=mean(exploration.guild.RA$mineral.soil.C.N), 
                     Distance=mean(exploration.guild.RA$Distance),
                     slope="no",
                     decomposing.carcass="no")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")

# saprotroph diversity increased with decomposing carcass and organic soil C:N, decreased with mineral soil C:N and distance

#######################
# SYMBIOTROPHS #
############################

# relative abundance
all.parms <- lm(symbiotroph ~ organic.soil.C.N + mineral.soil.C.N + organic.soil.GWC + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 +
                  slope + 
                  #ph + 
                  ln(Distance) + carcass.load +
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(formula = symbiotroph ~ paper.birch + moss + fern + horsetail + hogsweed + willow + azalea + ln(Distance) + carcass.deposition + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# symbiotroph relative abundance decreased with distance and carcass deposition

# richness
all.parms <- lm(symbiotroph.richness ~ organic.soil.C.N + mineral.soil.C.N + organic.soil.GWC + slope + ln(Distance) + carcass.load +
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 + ph +
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(formula = symbiotroph.richness ~ ln(Distance) + paper.birch + vaccinium + heather + eriophorum + horsetail + alder + hogsweed + azalea + carcass.deposition + 
                          (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)
# symbiotroph richness decreased with mineral soil C:N and carcass deposition, increased with organic soil C:N

# diversity 
all.parms <- lm(symbiotroph.diversity ~ organic.soil.C.N + mineral.soil.C.N + organic.soil.GWC  + slope + ln(Distance) + carcass.load + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 + ph + 
                  carcass.deposition + decomposing.carcass, 
                data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

# symbiotroph diversity decreased with mineral soil C:N
# without plants, decreased with carcass deposition and mineral soil C:N

#########################################
# SHORT DISTANCE TYPES # 
#########################################

# step model for relative abundance of short distance types
options(na.action = "na.fail")
all.parms <- lm(short.distance ~ organic.soil.C.N + mineral.soil.C.N + GWC + carcass.load + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 + ph + 
                  slope + ln(Distance) + carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
#dredge.model <- dredge(all.parms)
summary(step.model)
summary(dredge.model)
values$coefficients

step.model.lmer <- lmer(formula = short.distance ~ Dim1 + Dim2 + Dim3 + ln(Distance) + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# short distance relative abundance decreased with distance

all.parms <- lm(short.richness ~ organic.soil.C.N + mineral.soil.C.N + GWC + 
                  #carcass.load + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 + ph + 
                  slope + ln(Distance) + carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(formula = short.distance ~ mineral.soil.C.N + Dim4 + ph + ln(Distance) + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# short richness decreased with mineral soil C:N

# short diversity # 
exploration.guild.RA.short.div <- exploration.guild.RA %>% drop_na(short.diversity)
all.parms <- lm(short.diversity ~ organic.soil.C.N + mineral.soil.C.N + GWC + slope + ln(Distance) + carcass.load + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 + ph + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA.short.div)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(Distance=ln(mean(exploration.guild.RA$Distance)), 
                     mineral.soil.C.N=mean(exploration.guild.RA$mineral.soil.C.N),
                     Distance=mean(exploration.guild.RA$Distance),
                     carcass.deposition="yes",
                     decomposing.carcass="no")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")

step.model.lmer <- lmer(formula = short.distance ~ mineral.soil.C.N + ln(Distance) + 
                          carcass.deposition + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# short diversity decreased with mineral soil C:N and carcass deposition

#####################################
# MEDIUM DISTANCE FRINGE # 
#####################################

all.parms <- lm(medium.distance.fringe ~ organic.soil.C.N + mineral.soil.C.N + GWC + carcass.load +  
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 + ph + 
                  slope + ln(Distance) + carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(organic.soil.C.N=mean(exploration.guild.RA$organic.soil.C.N), 
                     mineral.soil.C.N=mean(exploration.guild.RA$mineral.soil.C.N), slope="yes",
                     carcass.deposition="yes")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")

step.model.lmer <- lmer(formula = medium.distance.fringe ~ organic.soil.C.N + mineral.soil.C.N + 
                          slope + carcass.deposition + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# medium distance fringe relative abundance increased with slope and organic soil C:N, and decreased with mineral soil C:N
# without plants, decreased with carcass deposition (p=0.07) and mineral soil C:N, and increased with slope and organic soil C:N

all.parms <- lm(medium.fringe.richness ~ organic.soil.C.N + mineral.soil.C.N + GWC + 
                  carcass.load + 
                  ph + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 +
                  slope + ln(Distance) + carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(organic.soil.C.N=mean(exploration.guild.RA$organic.soil.C.N), 
                     mineral.soil.C.N=mean(exploration.guild.RA$mineral.soil.C.N),
                     carcass.load=mean(exploration.guild.RA$carcass.load),
                     carcass.deposition="no")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")

step.model.lmer <- lmer(formula = medium.distance.fringe ~ organic.soil.C.N + carcass.load + 
                          ph + Dim1 + Dim2 + Dim3 + carcass.deposition + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# medium fringe richness increased with carcass deposition and distance

exploration.guild.RA.medium.div <- exploration.guild.RA %>% drop_na(medium.fringe.diversity)
all.parms <- lm(medium.fringe.diversity ~ organic.soil.C.N + mineral.soil.C.N + GWC + carcass.load + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 + ph + 
                  slope + ln(Distance) + carcass.deposition + decomposing.carcass, data=exploration.guild.RA.medium.div)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(Distance=mean(exploration.guild.RA$Distance),
                     decomposing.carcass="no")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")

step.model.lmer <- lmer(formula = medium.fringe.diversity ~ mineral.soil.C.N + carcass.load + 
                          Dim1 + Dim2 + Dim3 + ph + ln(Distance) + decomposing.carcass + (1|Stream:Transect), 
   data = exploration.guild.RA.medium.div)
summary(step.model.lmer)

  # medium fringe diversity increased with mineral soil C:N and decomposing carcass

##########################################################
# MEDIUM DISTANCE SMOOTH # 
##########################################################

exploration.guild.RA.hansen <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Hansen"),]
exploration.guild.RA.hansen.close <- exploration.guild.RA.hansen[which(exploration.guild.RA.hansen$Distance < 21),]

# relative abundance # 
all.parms <- lm(medium.distance.smooth ~ organic.soil.C.N + mineral.soil.C.N + GWC + carcass.load + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 +
                  slope + ph + 
                  ln(Distance) + carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(Distance=mean(exploration.guild.RA$Distance),
                     decomposing.carcass="no")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")

step.model.lmer <- lmer(formula = medium.distance.smooth ~ Dim1 + Dim2 + ln(Distance) + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# medium distance smooth relative abundance increase with decomposing carcass and decreased with distance

# richness # 
all.parms <- lm(medium.smooth.richness ~ organic.soil.C.N + mineral.soil.C.N + GWC + 
                  carcass.load + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 +
                  slope + ph + 
                  ln(Distance) + carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(Distance=ln(mean(exploration.guild.RA$Distance)), 
                     mineral.soil.C.N=mean(exploration.guild.RA$mineral.soil.C.N),
                     carcass.load=mean(exploration.guild.RA$carcass.load),
                     carcass.deposition="no")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")

step.model.lmer <- lmer(formula = medium.smooth.richness ~ mineral.soil.C.N + carcass.load + 
                          Dim1 + Dim2 + Dim3 + Dim4 + ph + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# medium smooth richness increased with carcass deposition (p=0.09)

# diversity #
exploration.guild.RA.medium.div <- exploration.guild.RA %>% drop_na(medium.smooth.diversity)
all.parms <- lm(medium.smooth.diversity ~ organic.soil.C.N + mineral.soil.C.N + GWC + slope + ln(Distance) + 
                  carcass.load +
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 + ph + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA.medium.div)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(GWC=mean(exploration.guild.RA$GWC),
                     mineral.soil.C.N=mean(exploration.guild.RA$mineral.soil.C.N),
                     decomposing.carcass="no")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")

step.model.lmer <- lmer(formula = medium.smooth.diversity ~ Dim1 + ph + decomposing.carcass + (1|Stream:Transect), data = exploration.guild.RA.medium.div)
summary(step.model.lmer)

# medium distance smooth diversity only affected by plants
# without plants, increased with decomposing carcass (p=0.07) and decreased with mineral soil C:n

##########################################################
# LONG DISTANCE # 
##########################################################

exploration.guild.RA.hansen <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Hansen"),]
exploration.guild.RA.hansen.close <- exploration.guild.RA.hansen[which(exploration.guild.RA.hansen$Distance < 21),]

# relative abundance # 
all.parms <- lm(long.distance ~ organic.soil.C.N + mineral.soil.C.N + carcass.load + 
                  mineral.soil.C.N + GWC + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 +
                  slope + ph + 
                  ln(Distance) + carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(Distance=mean(exploration.guild.RA$Distance),
                     decomposing.carcass="no")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")

step.model.lmer <- lmer(formula = long.distance ~ organic.soil.C.N + Dim3 + ph + decomposing.carcass + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# long distance decreased with distance and increased with decomposing carcass ()
# without plants, long distance increased with decomposing carcass

# richness # 
all.parms <- lm(long.richness ~ organic.soil.C.N + mineral.soil.C.N + GWC + carcass.load + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 +
                  slope + ph +
                  ln(Distance) + carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
newdata = data.frame(Distance=ln(mean(exploration.guild.RA$Distance)), 
                     mineral.soil.C.N=mean(exploration.guild.RA$mineral.soil.C.N),
                     carcass.deposition="no")
predict(step.model, newdata, se.fit=TRUE, interval="confidence")

step.model.lmer <- lmer(formula = long.richness ~ organic.soil.C.N + Dim1 + Dim3 + ph + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# long richness decreased with mineral soil C:N

# diversity # 
exploration.guild.RA.long.div <- exploration.guild.RA %>% drop_na(long.diversity)
all.parms <- lm(long.diversity ~ organic.soil.C.N + mineral.soil.C.N + GWC + carcass.load + 
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 +
                  slope + ph + 
                  ln(Distance) + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA.long.div)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(long.diversity ~ mineral.soil.C.N + Dim4 + ph + 
                          ln(Distance) + carcass.deposition + (1|Stream:Transect), data = exploration.guild.RA.long.div)
summary(step.model.lmer)

# long distance diversity changed with plants and decreased with organic soil C:N


###########################################################################################
# INDIVIDUAL GENERA # 
###########################################################################################
library(phyloseqCompanion)
library(ggplot2)
library(Rmisc)
library(ggsignif)
library(ggpubr)
library(dplyr)
library("gridExtra")
library(microbiomeMarker)
library(tidyr)

# set scientific notation off
options(scipen=999)

# set directory
setwd("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code")

# load phyloseq variables
load("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code/phyloseq final variables.RData")

# remove bog sample
physeq_nobog <- subset_samples(physeq, Transect != "Bog")

sample_data(physeq_nobog)$organic.soil.GWC <- c(1:95)
sample_data(physeq_nobog)$Dim1 <- c(1:95)
sample_data(physeq_nobog)$Dim2 <- c(1:95)
sample_data(physeq_nobog)$Dim3 <- c(1:95)
sample_data(physeq_nobog)$Dim4 <- c(1:95)
sample_data(physeq_nobog)$Dim5 <- c(1:95)
sample_data(physeq_nobog)$organic.soil.C.N <- c(1:95)
sample_data(physeq_nobog)$mineral.soil.C.N <- c(1:95)
sample_data(physeq_nobog)$carcass.deposition <- c(1:95)
sample_data(physeq_nobog)$decomposing.carcass <- c(1:95)
sample_data(physeq_nobog)$ph <- c(1:95)
sample_data(physeq_nobog)$slope <- c(1:95)
sample_data(physeq_nobog)$carcass.load <- c(1:95)
sample_data(physeq_nobog)$alder <- c(1:95)

# add covariates to phyloseq object
for (i in 1:95){
  for (j in 1:95){
    if (exploration.guild.RA$Stream[i] == sample_data(physeq_nobog)$Stream[j] & exploration.guild.RA$Distance[i] == sample_data(physeq_nobog)$Distance[j] & 
        exploration.guild.RA$Side[i] == sample_data(physeq_nobog)$Side[j] & exploration.guild.RA$Transect[i] == sample_data(physeq_nobog)$Transect[j]){
      sample_data(physeq_nobog)$Dim1[j] <- exploration.guild.RA$Dim1[i];
      sample_data(physeq_nobog)$Dim2[j] <- exploration.guild.RA$Dim2[i];
      sample_data(physeq_nobog)$Dim3[j] <- exploration.guild.RA$Dim3[i];
      sample_data(physeq_nobog)$Dim4[j] <- exploration.guild.RA$Dim4[i];
      sample_data(physeq_nobog)$Dim5[j] <- exploration.guild.RA$Dim5[i];
      #sample_data(physeq_nobog)$alder[j] <- exploration.guild.RA$alder[i];
      sample_data(physeq_nobog)$organic.soil.C.N[j] <- exploration.guild.RA$organic.soil.C.N[i];
      sample_data(physeq_nobog)$mineral.soil.C.N[j] <- exploration.guild.RA$mineral.soil.C.N[i];
      sample_data(physeq_nobog)$organic.soil.GWC[j] <- exploration.guild.RA$organic.soil.GWC[i];
      sample_data(physeq_nobog)$carcass.deposition[j] <- exploration.guild.RA$carcass.deposition[i];
      sample_data(physeq_nobog)$ph[j] <- exploration.guild.RA$ph[i];
      sample_data(physeq_nobog)$decomposing.carcass[j] <- exploration.guild.RA$decomposing.carcass[i];
      sample_data(physeq_nobog)$slope[j] <- exploration.guild.RA$slope[i];
      sample_data(physeq_nobog)$carcass.load[j] <- exploration.guild.RA$carcass.load[i];}
    }
  }
  

# read in soil data
soil.data <- read.csv("soil.data.final.csv")

# use relative abundance
physeq.prop <- transform_sample_counts(physeq_nobog, function(x) x/sum(x))

# agglomerate by genus
genus.glom <- tax_glom(physeq.prop, taxrank = 'Genus')

# or agglomerate by species
species.glom <- tax_glom(physeq.prop, taxrank = 'Species')

# get genus abundance matrix from phyloseq object
dat_genus <- psmelt(genus.glom)
dat_genus$Genus <- as.character(dat_genus$Genus)
dat_genus$Genus <- dat_genus$Genus %>% replace_na('Unclassified')
genus.abund <- subset(dat_genus, select = c("Sample", "Abundance", "Stream", "Transect", "Side", "Distance", 
    "organic.soil.C.N", "mineral.soil.C.N","organic.soil.GWC","carcass.deposition", "ph", "decomposing.carcass", 
    "Genus", "slope", "carcass.load", "Dim1", "Dim2", "Dim3", "Dim4", "Dim5"))

# create unique number for each genera
genus.abund <- genus.abund %>%
  group_by(Genus) %>%
  mutate(ID = cur_group_id())

genus.abund <- as.data.frame(genus.abund)

# long distance
Paxillus <- genus.abund[which(genus.abund$Genus=="Paxillus"),]
Boletus <- genus.abund[which(genus.abund$Genus=="Boletus"),]
Alpova <- genus.abund[which(genus.abund$Genus=="Alpova"),]
Xerocomus <- genus.abund[which(genus.abund$Genus=="Xerocomus"),]

# medium distance fringe
Cortinarius <- genus.abund[which(genus.abund$Genus=="Cortinarius"),]
Amphinema <- genus.abund[which(genus.abund$Genus=="Amphinema"),]
Piloderma <- genus.abund[which(genus.abund$Genus=="Piloderma"),]
Sistotrema <- genus.abund[which(genus.abund$Genus=="Sistotrema"),]
Tricholoma <- genus.abund[which(genus.abund$Genus=="Tricholoma"),]
Lyophyllum <- genus.abund[which(genus.abund$Genus=="Lyophyllum"),]


# medium distance smooth #
Lactarius <- genus.abund[which(genus.abund$Genus=="Lactarius"),]
Tomentella <- genus.abund[which(genus.abund$Genus=="Tomentella"),]
Entoloma <- genus.abund[which(genus.abund$Genus=="Entoloma"),]
Pseudotomentella <- genus.abund[which(genus.abund$Genus=="Pseudotomentella"),]
Hydnum <- genus.abund[which(genus.abund$Genus=="Hydnum"),]
Thelephora <- genus.abund[which(genus.abund$Genus=="Thelephora"),]
Tomentellopsis <- genus.abund[which(genus.abund$Genus=="Tomentellopsis"),]

# short distance # 
Russula <- genus.abund[which(genus.abund$Genus=="Russula"),]
Naucoria <- genus.abund[which(genus.abund$Genus=="Naucoria"),]
Inocybe <- genus.abund[which(genus.abund$Genus=="Inocybe"),]
Clavulina <- genus.abund[which(genus.abund$Genus=="Clavulina"),]
Laccaria <- genus.abund[which(genus.abund$Genus=="Laccaria"),]
Tylospora <- genus.abund[which(genus.abund$Genus=="Tylospora"),]
Hebeloma <- genus.abund[which(genus.abund$Genus=="Hebeloma"),]
Meliniomyces <- genus.abund[which(genus.abund$Genus=="Meliniomyces"),]
Tuber <- genus.abund[which(genus.abund$Genus=="Tuber"),]
Leotia <- genus.abund[which(genus.abund$Genus=="Leotia"),]
Sebacina <- genus.abund[which(genus.abund$Genus=="Sebacina"),]
Alnicola <- genus.abund[which(genus.abund$Genus=="Alnicola"),]
Hymenogaster <- genus.abund[which(genus.abund$Genus=="Hymenogaster"),]
Pulvinula <- genus.abund[which(genus.abund$Genus=="Pulvinula"),]
Endogone <- genus.abund[which(genus.abund$Genus=="Endogone"),]
Trichophaea <- genus.abund[which(genus.abund$Genus=="Trichophaea"),]
Wilcoxina <- genus.abund[which(genus.abund$Genus=="Wilcoxina"),]
Cenococcum <- genus.abund[which(genus.abund$Genus=="Cenococcum"),]
Genabea <- genus.abund[which(genus.abund$Genus=="Genabea"),]
Otidea <- genus.abund[which(genus.abund$Genus=="Otidea"),]
Geopora <- genus.abund[which(genus.abund$Genus=="Geopora"),]
Amanita <- genus.abund[which(genus.abund$Genus=="Amanita"),]
Hygrophorus <- genus.abund[which(genus.abund$Genus=="Hygrophorus"),]
Helvellosebacina <- genus.abund[which(genus.abund$Genus=="Helvellosebacina"),]


# saprotrophs
Mortierella <- genus.abund[which(genus.abund$Genus=="Mortierella"),]
Oidiodendron <- genus.abund[which(genus.abund$Genus=="Oidiodendron"),]

Paxillus.hansen <- Paxillus[which(Paxillus$Stream=="Hansen"),]
Boletus.hansen <- Boletus[which(Boletus$Stream=="Hansen"),]
Xerocomus.hansen <- Xerocomus[which(Xerocomus$Stream=="Hansen"),]
Alpova.hansen <- Alpova[which(Alpova$Stream=="Hansen"),]
Russula.hansen <- Russula[which(Russula$Stream=="Hansen"),]

# PAXILLUS #
all.parms <- lm(Abundance ~ organic.soil.C.N + mineral.soil.C.N + organic.soil.GWC + slope + ln(Distance) + carcass.load+
                  Dim1 + Dim2 + Dim3 + Dim4 + Dim5 + ph +
                  decomposing.carcass + carcass.deposition, 
                data=Xerocomus)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(Abundance ~ organic.soil.C.N + organic.soil.GWC + 
                          slope + decomposing.carcass + (1|Stream:Transect), data = Xerocomus)
summary(step.model.lmer)

# write results to table
tidy_lmfit <- tidy(step.model)
tidy_lmfit$Species <- c("Helvellosebacina")
#tidy_lmfit_final <- tidy_lmfit
tidy_lmfit_final <- rbind(tidy_lmfit_final, tidy_lmfit)
write.csv(tidy_lmfit_final, "tidy_lmfit.csv")

# further explore richness of species at carcass experiment site
unique(all.otu.abundance.final$Species[which(all.otu.abundance.final$Genus=="Cortinarius" & all.otu.abundance.final$Stream=="Hansen" 
  & all.otu.abundance.final$Side=="SL" & all.otu.abundance.final$Distance < 7 & all.otu.abundance.final$Abundance > 0)])

unique(all.otu.abundance.final$Species[which(all.otu.abundance.final$Genus=="Cortinarius" & all.otu.abundance.final$Stream !="Hansen" 
   & all.otu.abundance.final$Side=="SR" & all.otu.abundance.final$Distance < 7 & all.otu.abundance.final$Abundance > 0)])

unique(all.otu.abundance.final$Species[which(all.otu.abundance.final$Genus=="Amphinema" & all.otu.abundance.final$Stream=="Hansen" 
   & all.otu.abundance.final$Side=="SL" & all.otu.abundance.final$Distance < 7 & all.otu.abundance.final$Abundance > 0)])

unique(all.otu.abundance.final$Species[which(all.otu.abundance.final$Genus=="Amphinema" & all.otu.abundance.final$Stream=="Hansen" 
   & all.otu.abundance.final$Side=="SR" & all.otu.abundance.final$Distance < 7 & all.otu.abundance.final$Abundance > 0)])

unique(all.otu.abundance.final$Species[which(all.otu.abundance.final$Genus=="Lyophyllum" & all.otu.abundance.final$Stream=="Hansen" 
   & all.otu.abundance.final$Side=="SL" & all.otu.abundance.final$Distance < 7 & all.otu.abundance.final$Abundance > 0)])

unique(all.otu.abundance.final$Species[which(all.otu.abundance.final$Genus=="Lyophyllum" & all.otu.abundance.final$Stream=="Hansen" 
   & all.otu.abundance.final$Side=="SR" & all.otu.abundance.final$Distance < 7 & all.otu.abundance.final$Abundance > 0)])

Russula.hansen <- Russula[which(Russula$Stream=="Hansen"),]
Russula.hansen.left <- Russula.hansen[which(Russula.hansen$Side=="SL"),]
Russula.hansen.right <- Russula.hansen[which(Russula.hansen$Side=="SR"),]

tgc <- summarySE(Russula.hansen.left, measurevar="Abundance", groupvars=c("Side", "Distance"))
p1 <- ggplot(tgc, aes(x=Distance, y=Abundance)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Paxillus")+
  theme_classic() + ylab("Relative Abundance") + scale_x_reverse(breaks=c(1,6,10,20, 40, 60, 100)) +
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9))+ ylim(0,0.7)

tgc <- summarySE(Russula.hansen.right, measurevar="Abundance", groupvars=c("Side", "Distance"))
p2 <- ggplot(tgc, aes(x=Distance, y=Abundance)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Paxillus")+
  theme_classic() + ylab("Relative Abundance") + scale_x_continuous(breaks=c(1,6,10,20, 40, 60, 100)) +
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + ylim(0,0.7)

p <- grid.arrange(p1, p2, nrow=1, ncol=2)
